using CSV
using DataFrames
using Plots
using Statistics
using XLSX

# 读取 Excel 文件中的数据

filename = "sample_data/20241202采样.csv"  # 替换为您的文件名
data = DataFrame(CSV.File(filename))


temperatures = data[:, 1]
measurements = data[:, 2]
point_numbers = 1:length(measurements)

## 提取各列数据
# 高精度电压表导出数据
# p_n | time | temp | measure 
# ————+——————+——————+————————
#   1 |  t1  |  T1  |  data1
#   2 |  t2  |  T2  |  data2
#  ...| .... | .... |   ...
#   n |  tn  |  Tn  |  datan
# point_numbers = data[:, 1]
# timestamps = data[:, 2]
# temperatures = data[:, 3]
# measurements = data[:, 4]

# adc采样导出数据格式
# measurements = data[:, 1]
# point_numbers = 1:length(measurements)
# timestamps = point_numbers = 1:length(measurements)
# temperatures = point_numbers = 1:length(measurements)

# 滑动滤波（Moving Average Filter）
window_size = 50  # 滑动窗口大小
moving_avg_filtered = [mean(measurements[max(1, i-window_size+1):i]) for i in 1:length(measurements)]

# 均值滤波（Mean Filter）
mean_filtered = mean(measurements)
mean_filtered_array = fill(mean_filtered, length(measurements))

function kalman_filter_step(x,P,z)
    # function kalman_filter_step(A,B,H,Q,R,x,P)
    dt = 0.02  # 时间步长
    A = 1.0   # 状态转移矩阵
    B = 0.0   # 控制输入矩阵
    H = 1.0   # 观测矩阵
    Q = 0.01  # 过程噪声协方差
    # R = 4.272e-7   # 测量噪声协方差
    R = 1   # 测量噪声协方差
    
    x_hat = x
    x_hat_minus = A * x_hat + B * 0
    P_minus = A * P * A' + Q
    
    # 更新步骤
    K = P_minus * H' / (H * P_minus * H' + R)
    x_hat = x_hat_minus + K * (z - H * x_hat_minus)
    P = (1 - K * H) * P_minus
    
    return x_hat_minus,x_hat,P
end

# 卡尔曼滤波（Kalman Filter）
function kalman_filter(z)
    x_hat_minus = 2.63394896
    x_hat = 2.63394896  # 初始状态估计
    P = 1.0      # 初始误差协方差
    
    filtered_values = zeros(length(z))
    
    for k in 1:length(z)
        # 预测步骤
        x_hat_minus,x_hat,P = kalman_filter_step(x_hat,P,z[k])
        filtered_values[k] = x_hat
    end
    
    return filtered_values
end

# # 卡尔曼滤波（Kalman Filter）
# function kalman_filter(z)
#     dt = 0.02  # 时间步长
#     A = 1.0   # 状态转移矩阵
#     B = 0.0   # 控制输入矩阵
#     H = 1.0   # 观测矩阵
#     Q = 0.01  # 过程噪声协方差
#     # R = 4.272e-7   # 测量噪声协方差
#     R = 1   # 测量噪声协方差
    
#     x_hat = 2.63394896  # 初始状态估计
#     P = 1.0      # 初始误差协方差
    
#     filtered_values = zeros(length(z))
    
#     for k in 1:length(z)
#         # 预测步骤
#         x_hat_minus = A * x_hat + B * 0
#         P_minus = A * P * A' + Q
        
#         # 更新步骤
#         K = P_minus * H' / (H * P_minus * H' + R)
#         x_hat = x_hat_minus + K * (z[k] - H * x_hat_minus)
#         P = (1 - K * H) * P_minus
        
#         filtered_values[k] = x_hat
#     end
    
#     return filtered_values
# end

kalman_filtered = kalman_filter(measurements)

# 绘制原始数据和滤波后的数据
filtered_data = DataFrame(
    Point_Number=point_numbers,
    Timestamp=timestamps,
    Temperature=temperatures,
    Original_Measurement=measurements,
    Moving_Avg_Filtered=moving_avg_filtered,
    Mean_Filtered=mean_filtered_array,
    Kalman_Filtered=kalman_filtered
)


output_filename = "result/Soft_Filter/filtered_data.xlsx"

if isfile(output_filename)
    println("文件 $output_filename 已存在，将删除并重新保存。")
    rm(output_filename)  # 删除现有文件
else
    println("文件 $output_filename 不存在，将直接保存。")
end

XLSX.writetable(output_filename, filtered_data)

println("滤波后的数据已保存到 $output_filename")
output_filename = "result/Soft_Filter/filtered_data.csv"

if isfile(output_filename)
    println("文件 $output_filename 已存在，将删除并重新保存。")
    rm(output_filename)  # 删除现有文件
else
    println("文件 $output_filename 不存在，将直接保存。")
end

CSV.write(output_filename, filtered_data)


# 输出提示信息
println("数据已保存到 $output_filename")

start_point = 100
end_point = 4000

p = plot(
    point_numbers[start_point:end_point], measurements[start_point:end_point],
    label="Original Measurement", linewidth=0.5, color=:blue,
    legend=:topright
)

# 添加移动平均值
plot!(p, point_numbers[start_point:end_point], moving_avg_filtered[start_point:end_point], label="Moving Average", linewidth=0.5, color=:red)

# 添加卡尔曼滤波后的数据
plot!(p, point_numbers[start_point:end_point], kalman_filtered[start_point:end_point], label="Kalman Filtered", linewidth=0.5, color=:green)

# 添加零线
# hline!([0], linestyle=:dash, color=:black, label="Zero Line")

# 设置标签和标题
xlabel!("Time")
ylabel!("Data")
title!("Filter Comparison")

# 保存图像
savefig("result/Soft_Filter/Filter_Compare.svg")

# 打印保存信息
println("Filter_Compare.png saved.")

start_point = 100
end_point = 4000

moving_avg_residuals = moving_avg_filtered-measurements
kalman_residuals = kalman_filtered-measurements

p = plot(
    point_numbers[start_point:end_point], 
    moving_avg_residuals[start_point:end_point], 
    label="Moving Average", linewidth=2, color=:red,
    legend=:topright
)


# 添加卡尔曼滤波后的数据
plot!(p, point_numbers[start_point:end_point], kalman_residuals[start_point:end_point], label="Kalman Filtered", linewidth=2, color=:green)

# 添加零线
hline!([0], linestyle=:dash, color=:black, label="Zero Line")

# 设置标签和标题
xlabel!("Time")
ylabel!("residuals")
title!("Filter residuals Comparison")

# 保存图像
savefig("result/Soft_Filter/Filter_Residuals_Compare.png")

# 打印保存信息
println("Filter_Residuals_Compare.png saved.")